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The stability features of steady states of the spherically symmetric 
bJ!)' Einstein- Vlasov system are investigated numerically. We find support 

for the conjecture by Zeldovich and Novikov that the binding energy 
maximum along a steady state sequence signals the onset of instability, 
' a conjecture which we extend to and confirm for non-isotropic states. 

^ . The sign of the binding energy of a solution turns out to be relevant 

for its time evolution in general. We relate the stability properties to 
the question of universality in critical collapse and find that for Vlasov 
matter universality does not seem to hold. 



1 Introduction 



For any given dynamical system the existence and stability of steady states 
is essential both from a mathematics and from a physics point of view. In 
this paper we investigate these questions by numerical means for the spher- 
ically symmetric, asymptotically flat Einstein- Vlasov system. This system 
describes a self gravitating collisionless gas in the framework of general rela- 
tivity. Here the matter is thought of as a large ensemble of particles, which 
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is described by a density function on phase space, and the individual par- 
ticles move along geodesies. The precise formulation of this system will be 
given in the next section; for further information on the Einstein- Vlasov 
system we refer to J]. In astrophysics the system is used to model compact 
star clusters and galaxies. In this context the stability question was first 
studied by Zeldovich e. a. in the sixties |3H1EZI- These authors character- 
ize a steady state by its central redshift and binding energy and conjecture 
that the binding energy maximum along a steady state sequence signals the 
onset of instability. For isotropic steady states Ipser [El and Shapiro and 
Teukolsky [HS] found numerical support for this conjecture. In our investi- 
gation we find that this conjecture holds for non-isotropic steady states as 
well, where the density on phase space depends on the particle energy and 
angular momentum. Moreover, we use three different kinds of perturbations 
(for details see Section 3) while in the previous studies mentioned above only 
one type of perturbation was implemented. We also find that the sign of 
the binding energy is crucial for the evolution of the perturbed solutions. A 
positive value implies that the solution is bound in the sense that not all 
matter can disperse to infinity, and in this case the perturbation of a stable 
state seems to lead to periodic oscillations. In the case of a negative value 
of the binding energy we observe that the solution disperses to infinity. 

Our initial motivation for studying the stability of steady states was its 
role in critical collapse. This topic started with the work of Choptuik 
where he studied the Einstein-scalar-field system. He took a fixed initial 
profile for the scalar field and scaled it by an arbitrary constant factor. This 
gives rise to a family of initial data depending on a real parameter A. It 
turned out that there exists a critical parameter A^, such that for A < A^, 
the corresponding solutions disperses as predicted by the theoretical results 
of Christodoulou |9 for small data, while for A > A^ the corresponding 
solutions collapse and produce a black hole in accordance with ^01 ■ The 
surprising result was that the limit of the mass M{A) of the black hole 
tends to zero for A ^ A^, so that in such a one-parameter family there are 
black holes with arbitrarily small mass. Choptuik found that this fact is 
related to the existence of self-similar solutions of the Einstein-scalar-field 
system, in particular, the critical solution is self-similar and universal, i.e., 
independent of the initial profile which determines the one-parameter family. 
Later on Choptuik, Chmaj, and Bizoh performed a similar investigation for 
the Einstein- Yang-Mills system [Sj. Here both cases lim^_»^^ M{A) = 
and limA^A^.,A>A, M{A) > were found and called type II and type I 
respectively. In the latter case there is a mass gap in the M(74)-curve. The 
possibility of type I behavior in this system was related to the existence 
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of the Bartnik-McKinnon solutions |H1 03, which are static. Again, the 
conjecture is that the critical solution is universal. 

As opposed to the field theoretic matter models mentioned above the 
Vlasov model is phenomenological, and in contrast to fiuid models several 
global results have been obtained for Vlasov matter. For the spherically sym- 
metric, asymptotically flat case it was shown in that sufficiently small 
initial data launch global, geodesically complete solutions which disperse for 
large times. It is also known that there do exist initial data, necessarily 
large, which develop singularities 29^. The proof relies on the Penrose sin- 
gularity theorem. There are no general results on the behavior of large data 
solutions yet, except for the following: If data on a hypersurface of constant 
Schwarzschild time give rise to a solution which develops a singularity after 
a finite amount of Schwarzschild time, then the first singularity occurs at 
the center of symmetry |26| . An analogous result where Schwarzschild time 
is replaced by maximal slicing has also been proved [22]. In |2] further re- 
sults on the global behavior of solutions for large data can be found. The 
transition between dispersion and gravitational collapse was numerically in- 
vestigated by Rein, Rendall, and Schaeffer [2Z| , and it was found that there 
is a mass gap in the M{A) curve. This result was later confirmed by Olabar- 
rieta and Choptuik JH]- In addition, the latter authors reported evidence 
that the mass gap is due to the presence of static solutions, and they found 
support that the critical solution is universal. 

In the present investigation we address the role of steady states in critical 
phenomena for the Einstein- Vlasov system and the question of universality 
by explicitly exploiting the fact that for this system the existence of an 
abundance of steady states is well-established |25j . and that these steady 
states can easily be computed numerically. Computing a steady state fs we 
consider the family Afg of initial data. Within every family of steady states 
given by a specific dependence on particle energy and angular momentum we 
find unstable ones which act as critical solutions: If they are perturbed with 
A > 1 they collapse to a black hole, if they are perturbed with ^ < 1 they 
either disperse or oscillate in a neighborhood of the steady state, depending 
on the sign of the binding energy. Due to the abundance of possible such 
dependences on particle energy and angular momentum there cannot be a 
universal critical solution in spherically symmetric collapse for the Einstein- 
Vlasov system. 

The paper proceeds as follows: In the next section we formulate the 
Einstein- Vlasov system, first in general coordinates and then in coordinates 
adapted to spherical symmetry. In j27j Schwarzschild coordinates were used. 
Here we use maximal areal coordinates, as was done in ^Hl- This has the 
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advantage that regions of spacetime containing trapped surfaces can be cov- 
ered. In Section 3 we discuss the numerical scheme which we use; it is a 
particle-in-ceh scheme of the type used for kinetic models in plasma physics. 
It should be noted that for the analogous scheme used in [27 a rigorous 
convergence proof has been established in j^H], and we conjecture that this 
can be done for the present scheme as well. In Section 4 we investigate the 
stability properties of certain steady states and compare our findings with 
the earlier work mentioned above. Section 5 is devoted to critical phenom- 
ena and non-universality, and in a final section we discuss the reliability of 
our code. 

Our main motivation for this numerical analysis is that it may lead to 
conjectures on the behavior of solutions of the Einstein- Vlasov system which 
may eventually be proven rigorously. 

2 The Einstein- Vlasov system 

We first write down the Einstein- Vlasov system in general coordinates on 
the tangent bundle TM of the spacetime manifold M. Following standard 
practice we normalize the physical constants to one. The system then reads 
as follows: 



Here / is the number density of the particles on phase space, E^^ and G^" de- 
note the Christoffel symbols and the Einstein tensor obtained from the space- 
time metric gab, \g\ denotes its determinant, T"'' is the energy-momentum 
tensor generated by /, x"- are coordinates on M, {x°',p^) the corresponding 
coordinates on the tangent bundle TM, Latin indices run from to 3, and 



is the rest mass of a particle at the corresponding point in phase space. We 
assume that all particles have rest mass 1 and move forward in time so that 
the distribution function / lives on the mass shell 




m = \gabP P 



,a 6:1/2 
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We consider this system in the spherically symmetric, asymptotically flat 
case and write the metric in the following form: 

ds^ = -{a^ + a^/3^)dt^ + 2a^pdtdr + a^dr^ + [dO^ + sin^ 6 dcj)"^) . 

Here the metric coefficients a, /3, and a depend on t € M and r > 0, a and 
a are positive, and the polar angles 6 G [0, n] and (j) G [0, 27r] parameterize 
the unit sphere. The radial coordinate r is thus the area radius. Let be 
the second fundamental form and deflne 

ra 

By imposing the maximal gauge condition, which means that each hyper- 
surface of constant t has vanishing mean curvature, we obtain the following 
field equations: 

a,. = ^a^K^ +47rraV+ 77-(l - a^), (1) 
2 2r 

3 

Kr = K — 4:TTaj, (2) 

r 

at = 2aaK + {aP)r, (3) 

a^^ = ar (— --] + ^ (2r— + a"^ - l) + 47ra^a(S - 3p). (4) 
\ a r r^^ \ a J 



The Vlasov equation takes the form 

dtf +(—-/?) drf + f-^ - 2aKw + ^]d^f = 0, (5) 
\ ae / \ a ar-^e 



where 



e = e{r,w,L) = y^l + w'^ + L/r"^. (6) 



The variables w and L can be thought of as the momentum in the radial 
direction and the square of the angular momentum respectively, see [201 
more details. The matter quantities are defined by 



(7) 



P{t,r) = — / ef{t,r,w,L)dLdw, 
^ J-ooJo 

j{i^r) = — / wf{t,r,w,L)dLdw, (8) 
^ J-ooJo 

pco poo 2 _j_ T IJI 

S(t,r) = ^ / ^ fit,r,w,L)dLdw. (9) 

"l^ J-oo Jo £ 
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We impose the following boundary conditions which ensure asymptotic flat- 
ness and a regular center: 

a{t, 0) = a{t, oo) = a{t, oo) = I. (10) 

The equations (^-Q together with the boundary conditions (fTU]) constitute 
the Einstein- Vlasov system for a spherically symmetric, asymptotically flat 
spacetime in maximal areal coordinates. 

Let us consider some simple properties of this system, in particular such 
as will be relevant for its numerical simulation; a more careful mathematical 
analysis of the system will be performed elsewhere First we note that 
the phase space density / is constant along solutions of the characteristic 
system 

a{T,r)w 



a{T, r)e 



P{r,r), (11) 



a^fr, r)e „ , , , , a(T,r)L 

^ = - 2a(r, r)K(r, r)w + / ' ' , (12) 

a[T,r) a{T,r)r'^€ 

L = (13) 

of the Vlasov equation. If r i— > {R,W, L){T,t,r,w, L) =: Z{T,t,z) de- 
notes the solution of the characteristic system with (R, W, L)(t, t, r, w, L) = 
(r, L) then 

/(t, r, w, L) = MR, W, L)(0, t, r, w, L)), 

o 

with / = /(O, •) the initial datum for /. Due to our choice of coordinates 
w,L in momentum space the characteristic flow is not measure preserving. 
Indeed, if by D we denote differentiation along a characteristic of the Vlasov 
equation, then 

D{f dL dw dr) = -( + f3r + —] f dL dw dr; (14) 
\a^e r ) 

the factor on the right hand side is the (r, L)-divergence of the right hand 
side of the characteristic system. Notice that by Eqn. Q this factor equals 
-Din a. If A C M+ X M X M+ is measurable and A{t) := Z(i,0,^) = 
{(i?, VF,L)(t,0,r,'u;,L) | ir,w,L) G A} then 

, f{t,z)dz= [ f{z) dz, 
A{t) J A a{t,R{t,0,z)) 

while 

^ a{t,r)fit,z)dz= I aiO,r)Xz)dz. 



I 

JA{t) 
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In particular, the total number of particles, which, since all particles have 
rest mass one, equals the rest mass of the system, is conserved: 

poo fOO fOO 

47rM / / a{t,r)f{t,r,w,L)dLdwdr = Mq. 

Jo J -co Jo 

Next we notice that Eqn. can be rewritten as 

(r^K,)^ = -Anr^aj. (15) 
Using Eqn. the second order equation @ can be rewritten as 

— Or ] = ^-nr^aaip + S) + Qr^aan^ . (16) 
a ) ^ 

The Hawking mass m is given by 



We also introduce the quantity 



1 

2 V 'a- 



^ •= n 1 1 - -2 ) ' 



and note that by Eqn. [i can be written in the form 

/x(t, r) = {^T^pit, s) + y^{t, s)^ s^ds. 

Assuming that the matter is compactly supported initially and hence also 
for later times Eqn. 1)151) implies that K;(t, r) ~ r^^ for r large. Hence the 
limits as r tends to oo of m and p are equal so that the ADM mass M can 
be written as 

M = y" ^47rp(t, r) + ^K;^(t,r)^ r^dr. 

The ADM mass is conserved. 

It is of interest to find out when trapped surfaces form, which means 
that 2m I r > 1, and in view of (|17|) this condition can be written as 
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Note also that 2m/ r < 1/2 implies that 2/i/r < 1/2. From ^H] we obtain 
the following inequalities 



a a 



< 1, 



a a 



< 1. 



(18) 



Finally we notice that the second order equation for a in the form H16|l can 
be integrated: 



ar{t,r) 



a{t,r) 



[4TTaa{p + S) + GaaK^) s^ds. 



(19) 



Thus a is monotonically increasing outwards, and from (jTU]) and H18() it 
follows that 

< a < 1, < a < 1. 



3 The numerical method 



o o 



Let us consider an initial condition /= /(r, L), define additional variables 
u > and (f) E [0, vr] by the relations 

2 2,-^ J. T 22-2i, 

u = w H — 7^, w = u cos (p, L = r u sm cp 



and assume that / vanishes outside the set {r,u,(j)) € [Rq,Ri] x [?7o)f^i] x 
[<I>Oi*l'i]- We will approximate the solution using a particle method. For a 
thorough treatment of particle methods in the context of plasma physics we 
refer to |^. To initialize the particles we take integers N^, N^, and define 



Ar = , Au = — — , A(p = 



Nr 



4> 



ri = Ro+{i-l) Ar, uj = Uq + ( j - I ] Au, 4>,, = + { k - ^ ] A(f>, 







fijk = f{ri,Uj,(j)k) inrfAr 2TrUjAu smcpkAcj), 
r%k = n, w%k = cos^fc, Lijk = {nuj sin^fc)^ . 

At this point it should be noted that f^-^^ contains the phase space volume 
element, a fact which is convenient when computing the induced components 
of the energy momentum tensor, but which has to be observed when this 
quantity is time-stepped. At each point in phase space with coordinates 
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,0 



■ tjk^'^ijk^ -^ijk we imagine a particle with weight f'^-j^ which is smeared out 
in the radial direction, i.e., with respect to r it is represented by a hat 
function of width 2Ar. Once the particles are initialized the grid covering 
the support of the initial datum plays no further role. From these numerical 
particles, approximations are made of the quantities p, j, S defined in 0, 
(jHJ, (jU at the spatial grid points 



rj ■■= jAr, j 



0, 



note that this is now a new grid, which we extend from r = to the radius 
of the spatial support of the matter quantities. The field equations and 
H15() . which do not contain the metric coefficient a, can now be integrated 
on the support of the matter from the origin outward, using the boundary 
conditions UlUj) and (r'^K)|^=Q = 0; notice that k is then also known outside 
the support of the matter. 

To obtain an approximation of the metric coefficient a at the grid points 
rj we discretize Eqn. in the following way: We write 



and from what is already computed we obtain approximations for a, k, p, S 
both at the grid points rj and by interpolation at the intermediate points 
I'n^i, which we denote by 

■1^ y. 



aj, aj^i, 



Using the approximations 

„2 



-a, 



{rj 



and 



1 

' Ar 



a-i 

J 1 



a 



^, ar{r-i) 



a,- 



«7-l 



Ar 2 Ar 

we obtain the following linear system for the approximations aj for a{rj) 
on the grid points: For j = 0, . . . , N — 1, 

Uj+i H -aj-i - 



r 



7 + i 

2 _|_ -' 2 



a 



'3+h 



2 



+ [6^2 + 47r(pj + Sj)] r]{Ar)'^aj 



0, 
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r 



and the approximation 



for r large, in particular, r well outside the support of the matter. The latter 
approximation is motivated by the known representation of the asymptoti- 
cally flat Schwarzschild solution in maximal areal coordinates. The approxi- 
mation M* for ^(t, r) at large values of r can be computed from the already 
approximated quantities p and k. 

The linear system for aj is tridiagonal, obviously diagonally dominated, 
and it can easily be solved. Thus we now have approximations for all the 
field quantities on the spatial grid points rj. In passing we note that in ^H] 
the field equation Q was discretized directly resulting again in a tridiagonal 
system for a^-, but our approach of discretizing the rewritten version ()16() 
instead seems to perform better numerically. 

To perform the time step we propagate the numerical particles according 
to the characteristic system of the Vlasov equation (fTT|) . (|T2|) . (|T3|) . To do 
so we interpolate the field quantities to particle locations and use a sim- 
ple Euler time stepping method to define the new phase space coordinates 
''ijfc' ^ijfc' ^\jk °f numerical particle with label ijk. We still need to prop- 
agate the phase space volume element along the characteristics. Discretizing 
(|14j) in time we obtain the relation 



At 



1 



(/ ijk f ijk) 



a^'e r J 



which we use to compute jr., and one time step is complete. 



4 Stability issues for steady states 



For static solutions, /3 = k = j = 0so that 
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and the field equation decouples and is solved by 



„(,)=( l-^V". (20) 



The metric coefficient a is determined by the equation 



where a' = Or = da/dr and 

p{r) := ^ / — f{r,w,L)dLdw 
T J-ooJO e 

is the radial pressure. Eqn. (|21() is the rr-component of the Einstein equa- 
tions, and a lengthy computation using the Tolman-Oppenheimer-Volkov 
equation 

a' 2 

p' = {P + P) - -(p-Pt) 

a r 

shows that the second order field equation Q follows; 

PTir) ■■= ^ / -^f{r,w,L)dLdw 
■^f J-ooJo 

is the tangential pressure. It is easy to check that in addition to L also the 
particle energy 

E := a{r)^ 1 + w"^ + ^ = a(r)e 

is constant along characteristics in the static case. Hence the ansatz 

f{r,w,L) = ^{E,L), (22) 

satisfies the static Vlasov equation. By substituting it into the definitions for 
p and p these quantities become functionals of a, and the static Einstein- 
Vlasov system is reduced to the equation (|2T|) with (|20j) substituted in. 
Given some ansatz function $ and prescribing some value for a(0), the 
static system can therefore easily be solved numerically, by integrating (|2H) 
from r = outward, in particular if $ is such that the resulting dependence 
of p and p on a can be computed explicitly. This is the case for the power 
law 

f{r,w,L) = ^{E,L) = [E^-EfAL-L^)\. (23) 
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Here / > —1/2, k > —1, Lq > 0, Eq > is the cut-off energy, and x+ := 
max{x,0}. In the Newtonian case with I = Lq = this ansatz leads to 
steady states with a polytropic equation of state. Note that by taking Lq > 
there will be no matter in the region 



since there necessarily E > Eq and / vanishes. We call such configurations 
static shells. The existence of static shells with finite mass and finite ex- 
tension has been proved in '2T". It is numerically convenient to work with 
shells since potential difficulties in treating r = are avoided, and we will 
only consider shells here. This is also motivated by the fact that in the 
numerical experiments performed on critical phenomena in |^ the initial 
data for the matter had such a shell structure so that static shells are natu- 
ral candidates for the critical solutions. It should be pointed out that there 
are static solutions which do not have the form ()22[) . cf. [31j . 
In our simulations we have studied four cases: 



Case 


1: 


k = 


0, 


I 


= 0, 


Case 


2: 


k = 


0, 


I 


= 1/2, 


Case 


3: 


k = 


1, 


I 


= 1/2, 


Case 


4: 


k = 


0, 


I 


= 3/2. 



Having chosen k and / we then numerically construct static solutions to the 
Einstein- Vlasov system as indicated above, by specifying values on Eq, Lq 
and a(0). The resulting metric coefficient will a-priori not satisfy the bound- 
ary condition a{oo) = 1, but by shifting both a and £^0 appropriately a 
steady state which satisfies the boundary conditions is obtained. The dis- 
tribution function fg of the steady state is then multiplied by an amplitude 
A, so that a new, perturbed distribution function is obtained. This is then 
used as initial datum in our evolution code. Accordingly, if we choose A = 1 
then the initial datum is exactly the steady state, and a good test of our 
code is to check how much it deviates from being static in the evolution. 
We find that such initial data are tracked extremely well which is a very 
satisfying feature of our code. Of course, for unstable steady states the nu- 
merical errors introduced will make the solution drift off after some time 
which can be made longer by increasing the number of numerical particles 
and the number of time steps. 

For k and / fixed we characterize each steady state by its central red 
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shift Zc and its fractional binding energy £^6, which are defined by 



1, 



Mo' 



where ei, = Mq — M. 



The central redshift is the redshift of a photon emitted from the center and 
received at infinity, and the binding energy is the difference of the rest 
mass and the ADM mass and is a conserved quantity. In Figure 1 and 
Figure 2 below the relation between the fractional binding energy and the 
central redshift is given for Case 1 and Case 4 with Lq = 0.1. 



0.2 



-0.2 



-0.4 



0.2 0.4 0.6 0.6 



1 

Zc 



1.2 1.4 1.6 1.1 



Figure 1: Case 1 with Lq = 0.1 

The relevance of these concepts for the stability properties of steady 
states was first discussed by Zeldovich and Podurets [3H1 who argued that it 
should be possible to diagnose the stability from binding energy considera- 
tions. Zeldovich and Novikov j37j then conjectured that the binding energy 
maximum along a steady state sequence signals the onset of instability. Ipser 
jl7| and Shapiro and Teukolsky 1^ find numerical support for this conjec- 
ture and they also find that steady states with a central redshift above about 
0.5 result in collapse and are thus unstable. In both of these studies only 
isotropic steady states were considered, i.e. ^ = ^{E), whereas our study 
includes the dependence on L. Moreover, our algorithm is closely related to 
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0.1 




'o 0.5 1 1.5 2 2.5 3 3.5 4 4.5 

Zc 



Figure 2: Case 4 with Lq = 0.1 

the one in P7I for which convergence has been proved in [^. Of course, we 
also take advantage of the great progress of the computer capacity in that 
we can use a larger number of particles and considerably improve the reso- 
lution of the grid compared to what was possible in the earlier simulations. 
To get an indication of the accuracy of our simulation we compute at every 
time step the ADM mass M and the rest mass Mq, both of which should 
be conserved, and as long as no trapped surface has formed the errors are 
remarkably small. We will get back to a discussion of the numerical errors 
in the final section. 

Before describing the results of our simulations we mention that besides 
perturbations of a steady state fs of the form Afg we also considered per- 
turbations of the form fs{r-\-rsh,w, L) and fs{r,w + Wsh,L) where the state 
is shifted with respect to r by rgh or with respect to w by Wgh- It turns out 
that perturbations with yl > 1 or r^/j < or Wsh < 0, which in principle 
push the state towards collapse, and on the other hand perturbations with 
j4 < 1 or r^/j > or Wsh > 0, which in principle push the state towards 
dispersion, lead to the same qualitative features of the perturbed solution. 
In the following discussion we restrict ourselves to perturbations of the form 
Afs. 
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The general picture that arises from our simulations is summarized in 
Tables 1-4 which correspond to the four different cases we consider. The 
parameter Lq is the same in all cases, i.e. Lq = 0.1. 





Eb 


A<1 


A > 1 


0.24 


0.036 


stable 


stable 


0.30 


0.039 


stable 


stable 


0.39 


0.041 


stable 


stable 


0.43 


0.041 


stable 


unstable 


0.47 


0.040 


stable 


unstable 


0.52 


0.038 


stable 


unstable 


0.65 


0.027 


stable 


unstable 


0.82 


0.004 


stable 


unstable 


0.95 


-0.024 


unstable 


unstable 


1.1 


-0.070 


unstable 


unstable 


Table 1: Case 


1: A; = and / = 0. 


Zc 


Eb 


A<1 


A > 1 


0.21 


0.032 


stable 


stable 


0.34 


0.040 


stable 


stable 


0.39 


0.040 


stable 


stable 


0.42 


0.041 


stable 


unstable 


0.46 


0.040 


stable 


unstable 


0.56 


0.036 


stable 


unstable 


0.65 


0.029 


stable 


unstable 


0.82 


0.008 


stable 


unstable 


0.95 


-0.015 


unstable 


unstable 


1.20 


-0.078 


unstable 


unstable 



Table 2: Case 2: A; = and / = 1/2. 



If we first consider perturbations with j4 > 1 we find that steady states 
with small values on Zc (less than approximately 0.40 depending on which 
case we consider) are stable, i.e., the perturbed solutions stay in a neigh- 
borhood of the static solution as depicted in Figure 3. A more careful 
investigation of these perturbed solutions seems to indicate that they oscil- 
late in an (almost) periodic way, and we come back to this issue at the end 
of this section. For larger values of Zc the evolution leads to the formation 
of trapped surfaces and by the result in to the collapse to black holes, 
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Eh 


.4 < 1 


.4 > i 


U.Zi 


U.UoU 


stable 


stable 


0.33 


0.037 


stable 


stable 


0.38 


0.038 


stable 


stable 


0.42 


0.039 


stable 


stable 


0.045 


0.039 


stable 


unstable 


0.052 


0.037 


stable 


unstable 


0.77 


0.020 


stable 


unstable 


0.90 


0.004 


stable 


unstable 


1.0 


-0.01 


unstable 


unstable 


1.13 


-0.037 


unstable 


unstable 


Table 3: Case 3: A; = 1 and Z = 1/2. 






A< 1 


A > 1 


0.21 


0.032 


stable 


stable 


0.34 


0.039 


stable 


stable 


0.41 


0.040 


stable 


stable 


0.44 


0.040 


stable 


stable 


0.52 


0.038 


stable 


unstable 


0.69 


0.025 


stable 


unstable 


0.80 


0.014 


stable 


unstable 


0.92 


-0.004 


unstable 


unstable 


1.1 


-0.045 


unstable 


unstable 



Table 4: Case 4: k = and I = 3/2. 



as depicted in Figure 4. As a matter of fact we do check if null geodesies 
that start after a certain time from the center are caught so that they can- 
not escape to infinity, and our results always support that there is an event 
horizon. 

Hence, for perturbations with A > 1 the value of Zc alone seems to 
determine the stability features of the steady states. By plotting the curve. 
Eh versus Zc, for a shorter interval of Zc to get better resolution (Figures 5 
and 6 below) we find that the conjecture by Novikov and Zeldovich, that 
the maximum of along a sequence of steady states signals the onset of 
instability, is true in a numerical sense also for states that depend on the 
angular momentum L. Having a second look at Figure 2 one might think 
that along that curve another change of stability behavior might occur at 
Zc ~ 2, but we found no indication of this. The reason for the qualitative 
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t=2T/8 




Figure 4: = 1.14, Eb = -0.08, A = 1.01, T = 10.0 
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Zc 



Figure 5: Case 1 with Lq = 0.1 



0.06 




Zc 



Figure 6: Case 4 with Lq = 0.1 
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difference between Figure 1 and Figure 2 and its possible consequences will 
be investigated elsewhere jl]. 

The situation is quite different for perturbations with A < 1. The crucial 
quantity in this case is the fractional binding energy Ef,. Consider a steady 
state with Ef, > and a perturbation with A < 1 but close to 1 so that 
the fractional binding energy remains positive. Then the perturbed solution 
drifts outwards and then turns back and reimplodes and comes close to its 
initial state, and then continues to expand and reimplode and thus oscillates. 
This is depicted in Figure 7, see also the end of this section. 



t=0 



t=T/8 



t=2T/8 




Figure 7: = 0.47, Eh = 0.04, A = 0.99, T = 90.0 

In I^S] it is stated that if Ef, > the solution must ultimately reimplode. 
We are not aware of any precise mathematical formulation (or proof) of such 
a statement but our simulations support that it is true. For negative values 
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of Eh the solutions with A <1 disperse to infinity as depicted in Figure 8. 



o. 1 
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t=2T/8 




t=3T/8 



2 4 
t=4T/8 



2 4 
t=5T/8 



A 



2 4 

t=6T/8 



2 4 

t=7T/8 



2 4 
t=T 




Figure 8: Zc = 1.14, = -0.08, A = 0.99, T = 14.0 

A simple argument which relates the binding energy to the question 
whether a solution disperses or not at least for the case where the spatial 
support is a shell is the following: Consider a solution which has an expand- 
ing vacuum region of radius R{t) at the center with R{t) — oo for t ^ oo. 
Such a solution disperses in a strong sense, and we claim that this implies 
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that Mq < M, i.e. £^5 < 0. To see this, observe that 

f'OQ f'OO f'OO 

Mo = 47r2/ / / a{t,r) f{t,r,w,L)dLdw dr 
jRit) J -00 Jo 



Att^ I "^'^''' j j j f{t,r,w,L)dLdwdr 



R{t) J-00 JO 



-1/2 



< 1 



1--— 4^2/ / / ef{t,r,w,L)dLdwdr 

V ^{i)/ JR{t)J-ooJO 

so that with t ^ oo necessarily Mq < M as claimed. 

For the Vlasov-Poisson system, which is the Newtonian limit of the 
Einstein- Vlasov system [23] a rigorous stability theory for steady states of 
the form analogous to has been developed in recent years, based on 
variational techniques, cf. ^1 El El El 1121 • "^^^ result is that the stabil- 
ity properties of the steady state are essentially determined by the ansatz 
function $ in (|22() . in particular, all the polytropic steady states resulting 
from an ansatz of the form ((221) with < k < 7/2, 1 = 0, Lq = are non- 
linearly stable. The above discussion shows that an analogous result does 
not hold for the Einstein- Vlasov system, since there the stability depends 
in addition on the central redshift and the binding energy, even if k and I 
are fixed. We note at this point that the stability of shell like steady states 
has not yet been systematically investigated for the Vlasov-Poisson system. 
But we would argue that this does not affect the above statement, since on 
the one hand it is reasonable to expect that the stability results for shells 
in the Newtonian case will be analogous to the case Lq = 0, an expectation 
supported by recent results in this direction [32], and on the other hand 
we do not expect our numerical stability results for the relativistic case to 
depend on the shell structure of the steady states. A clear indication that 
stability results based on variational techniques do not carry over to the 
relativistic case without making additional restrictions are the first results 
in this direction in |36j . 

Our investigation indicates that there are (at least) three qualitatively 
different types of behavior which can result from the perturbation of a steady 
state. Perturbing an unstable state with Zc > 0.4 using A > 1 leads to 
gravitational collapse and a black hole. Perturbing an unstable state with 
Ej, < using A < 1 leads to dispersion. In the other, stable cases Zc < 0.4 
and ^ > 1 or > and ^ < 1 the perturbation leads to an oscillatory 
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behavior. Whether the perturbed solution is indeed time-periodic is hard 
to decide numerically. One way to shed some light on this question is to 
plot a{t,0) or the inner and outer radius rjnin(t) and r^ax{t) of the matter 
support. The results are shown in Figures IHl and ^1 for a steady state with 
k = 0,1 = 1/2 (Case 2) and Lq = 0.1, = 0.21, Eh = 0.03, perturbed 
with A = 1.02. The oscillatory behavior with a clear time period is quite 

0.86 I , , , , , 1 



0.85 - 




0.79 I • • • • • 1 

50 100 150 200 250 300 

t 



Figure 9: = 0.21, Eb = 0.03, A = 1.02, T = 300.0 

striking. We found that all the solutions arising by a small perturbation of 
a fixed steady state seem to have roughly the same period, independently of 
A. We do at this point venture no interpretation of the slow upward drift 
of a{t,0) or of rmax(i)- 

5 Critical phenomena and non-universality 

Critical collapse for the Einstein- Vlasov system has previously been studied 
numerically in |27 1 ll9 | l35]. By specifying an initial datum and then varying 
its amplitude A, a critical value Ac of A is found in the sense that if A > Ac 
the evolution of the initial datum will form a trapped surface and collapse 
to a black hole, and if A < Ac no black hole will ever form. In the previous 
studies referred to above it was stated that the solutions in the subcritical 
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Figure 10: = 0.21, E^, = 0.03, A = 1.02, T = 300.0 



case (i.e. A < Ac) will disperse to infinity. However, as we saw in the 
previous section this is not quite true since there are (at least) two possible 
scenarios: Either the solution will disperse to infinity or it will be bound 
depending on the sign of the binding energy e^. Nevertheless, the value 
Ac separates two distinct behaviors, as was first discovered by Choptuik 
who studied the Einstein-scalar field system. From his work a new topic 
developed which goes under the name of critical phenomena, cf. [EI for a 
review. Two types of matter are distinguished in this context. Vlasov matter 
is said to be of type I whereas a scalar field and a perfect fluid are of type II. 
One fundamental difference between matter of type I and type II is that the 
graph which describes the mass of the black hole (which of course is zero if 
A < Ac) as a function of A is discontinuous at A = Ac for type I matter, and 
continuous for type II matter. For type I matter there is thus a mass gap 
and no possibility to get a black hole with arbitrarily small mass as in the 
case of type II matter. The characteristic behavior of the critical solution 
itself has been carefully investigated for type II matter. A lot of evidence 
has been found that it is self-similar and universal, where the latter property 
means that it is independent of the initial datum, cf. ^2j. For type I matter 
the general conjecture is that the critical solution is static (and unstable) 
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and universal, see (71. This conjecture has been investigated in 19] and [SS] 
for Vlasov matter, and evidence that the critical solution is static was found. 
In these investigations it is also claimed that some support for universality is 
obtained, but it is concluded that further studies need to be done. From the 
previous section it is clear that only the static property, and not universality, 
can be a genuine feature of the critical solution for Vlasov matter. Indeed, 
by picking /o = /s, where fs is the distribution function of an unstable 
static solution we find Ac = 1 as the critical amplitude and fs as the critical 
solution. Since there are infinitely many unstable static solutions, and any 
one will do as critical solution, universality is contradicted. More precisely, 
within any family of steady states (which we have investigated), specified 
by some choice of k and / in (|23() . there are infinitely many unstable ones 
each of which acts as a critical solution. And the class of steady states of 
the form H23() is only a small subset of all possible forms of steady states 
obtainable by the ansatz (|22|) , cf . [2Sj . 

If we turn to the other conjectured property, that the critical solution 
is static, our results strongly support that this is the case. We start from 
subcritical initial data and tune the amplitude to a very high degree of 
precision to get Ac- We then find that in the evolution the metric coefficients 
become more or less frozen after some time and stay so during a considerably 
long time period before the solution starts to drift off outwards. The length 
of the latter time period depends on how carefully we have calibrated Ac- 
By a considerably long time period we mean considerable in comparison to 
a typical dynamical time scale, e. g. a typical time to form a trapped surface 
or the typical period time of an oscillation for a bound state. This is further 
discussed in the next section on the error estimates. 

6 The numerical errors 

To get an indication of the validity of our simulations we compute at every 
time step the ADM mass M{t) and the rest mass Mo(t), which are conserved 
along true solutions of the system. We define the errors e[M] and e[Mo] by 



where M = M(0) and Mq = Mo(0) are given by the initial datum. As we 
have explained in previous sections, the evolution of an initial datum will 
either form a trapped surface and collapse, will be bound and never collapse, 
or it will disperse to infinity. In the latter two situations the errors in our 



e[M]{t) 



\M{t) - M 
M 



e[Mo]{t) 



\Mo{t)-Mo\ 
Mo 
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simulations are remarkably small as can be seen in the table below. On the 
other hand when a trapped surface has formed a black hole singularity will 
develop (cf. The coordinates that we are using are believed to avoid 

the singularity so that solutions will be regular for all times. However, the 
solutions become more and more peaked at a certain more or less frozen 
radius. Of course a very fine grid is needed to resolve such sharp peaks. 
We are using a fixed grid which does not change with time, and since the 
solutions get more and more peaked as time goes on our grid is not sufficient 
after a certain time of the collapse. This results in an essential increase of 
the errors as seen in Table 5 in the case of collapsing initial data. 







A 


TTS 


PT 


TF 


e\M\{TF) 


e[Mo](TF) 


0.30 


0.039 


1.01 






50 


0.00035 


0.000036 


0.47 


0.039 


0.99 




60 


150 


0.0013 


0.000023 


1.14 


-0.080 


0.99 






20 


0.00015 


0.00011 


1.14 


-0.080 


1.01 


3.6 




10 


0.080 


0.040 



Table 5: Numerical errors 



For each of the simulations in Table 5 about 14 000 particles have been 
used. This gives convenient running times of a few minutes, but one can of 
course use a much larger number of particles and for certain runnings we 
have used up to 10^ particles which gives running times of several hours on 
a reasonably modern PC. To get an idea of the length of the final running 
time, TF, we have in the relevant cases included the time when a trapped 
surface forms, TTS, and the period time, PT, of an oscillation. In the 
collapsing situation (row 4) where the errors are of a different magnitude we 
also computed the error at t = 2 * TTS = 7.2 and found e[M](7.2) = 0.024 
and e[Mo](7.2) = 0.046, in order to get an indication of the growth of the 
errors. In conclusion, Table 5 clearly demonstrates that we can keep M and 
Mq conserved to a very high precision in the evolution as long as we are not 
considering collapsing solutions. In the latter case the errors are still very 
reasonable for a considerable running time compared to the time when a 
trapped surface forms. 



References 



[1] H. Andreasson, The Einstein- Vlasov system/kinetic theory. Liv- 
ing Rev. Relativ. 8 (2005). 



26 



[2] H. Andreasson, On global existence of the spherically symmet- 
ric Einstein- Vlasov system in Schwarzschild coordinates. Indiana 
University Math. J. To appear. 

[3] H. Andreasson, G. Rein, The spherically symmetric Einstein- 
Vlasov system in maximal areal coordinates. In preparation. 

[4] H. Andreasson, G. Rein, On the steady states of the spherically 
symmetric Einstein- Vlasov system. In preparation. 

[5] R. Bartnik, J. McKiNNON, Particlclikc solutions of the Einstein- 
Yang-Mills equations. Phys. Rev. Lett. 61, 141-143 (1988). 

[6] C. K. BiRDSALL, A. B. Langdon, Plasma Physics via Computer 
Simulation. New York: McGraw-Hill 1985 

[7] M. W. Choptuik, Universality and scaling in the gravitational 
collapse of a scalar field, Phys. Rev. Lett. 70, 9-12 (1993). 

[8] M. W. Choptuik, T. Chmaj, P. Bizon, Critical behaviour in 
gravitational collapse of a Yang-Mills field, Phys. Rev. Lett. 77, 
424-427 (1996). 

[9] D. Christodoulou, The problem of a sclf-gravitating scalar field. 
Commun. Math. Phys. 105, 337-361 (1986). 

[10] D. Christodoulou, A mathematical theory of gravitational col- 
lapse. Commun. Math. Phys. 109, 613-647 (1987). 

[11] M. Dafermos, a. D. Rendall, An extension principle for 
the Einstein- Vlasov system in spherical symmetry, Ann. Henri 
Poincare 6, 1137-1155 (2005) 

[12] C. GUNDLACH, Critical phenomena in gravitational collapse, Liv- 
ing Rev. Relativ. 1999-4, (1999). 

[13] Y. Guo, G. Rein, Stable steady states in stellar dynamics. Arch. 
Rational Mech. Anal. 147, 225-243 (1999) 

[14] Y. Guo, G. Rein, Existence and stability of Camm type steady 
states in galactic dynamics. Indiana University Math. J. 48, 1237- 
1255 (1999) 

[15] Y. Guo, G. Rein, Isotropic steady states in galactic dynamics. 
Commun. Math. Phys. 219, 607-629 (2001) 



27 



[16] Y. Guo, G. Rein, Stable models of elliptical galaxies. Monthly 
Not. Royal Astr. Soc. 344, 1396-1406 (2003) 



[17] J. R. Ipser, Relativistic, spherically symmetric star clusters. III. 
Stability of compact isotropic models. The Astrophys. J. 158, 17- 
43 (1969) 

[18] E. Malec, N. 6 MURCHADHA, Optical scalars and singularity 
avoidance in spherical spacetimes, Phys. Rev. D. 50, 6033-6036 
(1994). 

[19] I. Olabarrieta, M. W. Choptuik, Critical phenomena at the 
threshold of black hole formation for collisionless matter in spher- 
ical symmetry, Phys. Rev. D. 65, 024007 (2002). 

[20] G. Rein, The Vlasov-Einstein system with surface symmetry, Ha- 
bilitationsschrift, Munich, (1995). 

[21] G. Rein, Static shells for the Vlasov-Poisson and Vlasov-Einstein 
systems. Indiana University Math. J. 48, 335-346 (1999). 

[22] G. Rein, Nonlinear stability of Newtonian galaxies and stars from 
a mathematical perspective. In: Nonlinear Dynamics in Astron- 
omy and Physics, Annals of the New York Academy of Sciences 
1045, 103-119 (2005) 

[23] G. Rein, A. D. Rend all, Global existence of solutions of the 
spherically symmetric Vlasov-Einstein system with small initial 
data, Commun. Math. Phys. 150, 561-583 (1992). Erratum: Com- 
mun. Math. Phys. 176, 475-478 (1996). 

[24] G. Rein, A. D. Rend all. The Newtonian limit of the spherically 
symmetric Vlasov-Einstein system, Commun. Math. Phys. 150, 
585-591 (1992). 

[25] G. Rein, A. D. Rendall, Compact support of spherically sym- 
metric equilibria in non-relativistic and relativistic galactic dynam- 
ics. Math. Proc. Camb. Phil. Soc. 128, 363-380 (2000) 

[26] G. Rein, A. D. Rendall, J. Schaeffer, A regularity theorem 
for solutions of the spherically symmetric Vlasov-Einstein system, 
Commun. Math. Phys. 168, 467-478 (1995). 



28 



[27] G. Rein, A. D. Kendall, J. Schaeffer, Critical collapse of 
collisionless matter — a numerical investigation, Phys. Rev. D 58, 
044007 (1998). 

[28] G. Rein, T. Rodewis, Convergence of a particle-in-cell scheme 
for the spherically symmetric Vlasov-Einstein system, Indiana Uni- 
versity Math. J. 52, 821-862 (2003). 

[29] A. D. Rend ALL, Cosmic censorship and the Vlasov equation. 
Class. Quantum Grav. 9, L99-L104 (1992). 

[30] A. D. Rend ALL, Cosmic censorship and the Vlasov equation. 
Class. Quantum Grav. 9, L99-L104 (1992). 

[31] J. Schaeffer, A class of counterexamples to Jeans' theorem for 
the Vlasov-Einstein system. Comm. Math. Phys. 204, 313-327 
(1999). 

[32] A. SCHULZE, Nonlinear stability of stationary shells for the Vlasov- 
Poisson system. In preparation. 

[33] S. L. Shapiro, S. A. Teukolsky, Relativistic stellar dynamics 
on the computer. II — Physical Applications, Astrophysical J. 298, 
58-79 (1985) 

[34] J. A. Smoller, a. G. Wasserman, S.-T. Yau, J. B. McLeod, 
Smooth static solutions of the Einstein- Yang- Mills equations. Com- 
mun. Math. Phys. 143, 115-147 (1991). 

[35] R. Stevenson, The Spherically Symmetric Collapse of Collision- 
less Matter, M. Sc. Thesis, The University of British Columbia, 
(September 1 2005, DRAFT). 

[36] G. WOLANSKY, Static solutions of the Vlasov-Einstein system. 
Arch. Rational Mech. Anal. 156 205-230 (2001) 

[37] Ya. B. ZelDovich, I. D. NoviKOV, Relativistic Astrophysics 
Vol. 1, Chicago: Chicago University Press (1971) 

[38] Ya. B. ZelDovich, M. A. Podurets, The evolution of a system 
of gravitationally interacting point masses, Soviet Astronomy-AJ 
9, 742-749 (1965), translated from Astronomicheskii Zhumal 42, 
963-973 (1965) 



29 



